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This thesis attempts to develop a program for the Whirlwind 
eomputer operating in the interpretive mode for the flux distribution 
in a finite cylindrical reactor© 

The reactor chosen for investigation was a hypothetical organic 
moderated, highly enriched, cylindrical reactor having a power output 
capable of fulfilling one half the propulsive requirements of a 
nuclear powered naval cruiser© 

Heat transfer and basic physics calculations were accomplished 
to determine preliminary size and material requirement of such a 

reactor© 

The flux equations were derived using numerical integration 
methods suitable for computer programming© 
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Conclusion 



1) Tlie cod© as presently programmed is incorrect* 

2) The grid spacing as used is much too coarse* This 
difficulty is not inherent in the coded program* 

3) The mode of computer operation is unsuitable for an 
iterative process such as this because it results in 

• a code that progresses too slowly* 

Recommendations 

1) Prior to the use of any computer a formal or semi-fonnal 
course in programming should be congsleted* 

2) A grid spacing comparable to the mean free path of the 
neutron should be selected* 

3) The fastest possible mode of calculation on the computer 
available should be used* 

!i) For a coded program of this difficulty a faster machine 
with a much larger fast memory is very desirable* 
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INTR0I7JCTI0N 



The desirability of a vmiforai power distribution (P /P ) * 1 

max av 

has long been realized since it will not only decrease the mass of 
the fuel required for the attainment of criticality but also will 
increase the total power output by allowing the whole core instead 
of only a small fraction of the reactor to operate at the highest allow- 
able temperature* Operation at the highest possible temperature allowed 
by metallurgical considerations over the entire reactor will increase 
the average core life since more uniform bum up is attained* 

It is possible to achieve this uniform power distribution by 
various means* One of the first methods to be used was that of 
poisons* In this method long lasting poisons are placed in the 
reactor structure in the anticipated locations of high flux* By 
placing poison in these regions the peak flux is reduced to such a 
level that the power density is approximately uniform* This pro- 
oefiareais a very wasteful one from the point of view of neutron economy* 
T^e regions of the reactor that are contributing greatly to the overall 
neutron efficiency of the reactor sire penalized to the extent that 
they can be no more efficient than the most inefficient region* For 
power producing reactors, particularly mobile ones, this waste of 
neutrons is very undesirable* 

The use of poison is the first of two methods that ejtpteit a 
nonuniform macroscopic cross section to achieve the uniform 
distribution of power* The second consists in the use of a nonunj.form 
macroscopic fissiqn cross section* This method allows each regiop of 
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the reactor to operate most efficiently* The efficiency is determined 
by the position of the region within the reactor* The leakage from 
a reactor causes a reduction in the value of the flux at large 
distances from the center of the core* To compensate for this reduc- 
tion in flux the macroscopic fission cross section is increased in 
the i?egions far from the center and decreased in the central portions* 

To this end^ the fuel concentration is increased near the edges* 

Neutron economy can also be improved by surrounding the core with 
neutron reflector material* A light substance having a relatively 
low absorption cross section performs this function very well^ i*e*, 
heavy water or beryllium* 

Since the reactor to be investigated is intended for a hypothetical 
naval cruiser^ the use of poisons to achieve a uniform power distribu- 
tion is not desii’able^ for a compact, reactor is required* Thou^ the 
reactor must be compact it is desirable to install a reasonably thick 
reflector for the purpose of improving the neutron economy^ attenuating 
-thameutrons leaking out, and cooling of the reactor shell and 
pressure vessel* 

It would be possible to have a nonuniform macroscopic cross 
section of the reactor in each of its directions* For the spherical 
reactor the variation would be radial only, but the application of a 
spherical reactor to shipboard use is impractical* The norii^al con- 
figuration of spaces aboard a ship is rectangular and a spherical 
reactor does not adapt very well to this arrangement* However, a 
rectangular reactor allows efficient utilization of space aboard a 
ship, but from, the nuclear point of view this shape is not too desir- 
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able because leakage depends on the relationship of surface area to 
•volume. It can be seen that the neutron economy decreases in 
progression from the sphere to the cylinder to the rectangular 
parallepipedo 

With these considerations in mind a right cylindrical reactor 
shape vd.ll be investigated. With this shape it •would be possible 
to •vary the properties in tvo directions, radially and axially* 

While this goal is theoretically interesting, it is prohibi'ted 
when cost is considered* The reactor -that is tdbe investigated vd.ll have a 
urrLfom composi^tion in the axial direction -with fini-te radial segmen^bs 
differing in properties* 

For the shipboard use of a nuclear reactor power plant it is 
desirable •to have a plant as light and safe as can possibly be 
designed •within the space allovfed but still vri.thout incorporating 
undue quantities of the "exotic” and expensi^ve materials or coolan-ts* 

One possible solu-tion is the use of a highly enriched reactor using 
a coolant that permits hi^ operating temperatures at low pressures 
and does not become dangerovisly radioactive while passing thm 
the reac^bor core* 

A coolant that has a low •vapor pressure at high operating 
•temperature vd.ll allow a system pressure which ia much lower thah that 
eii5)loyed in the present day water cooled and moderated shipboard 
reac^bor installation* The use of such a coolant will not only result 
in much better steam condi^tions, but it vjill also make possible an 
overall weight reduc^tion, since a low pressure system can be used* 

The saving in weight is particularly significant in the pressure 
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If the coolant does not become radioactive in\4-ts passage thru 
the reactor core^ the shielding can be concentrated primarily around 
the reactor proper^ provided that impttrities in the coolant are 
kept at a very low concentration. Even with a small amoTxnt of 
irapuritieSg the radiation emitted by the coolant would be low 
enough to permit reduced thicknesses of shielding aroTind the prlraary 
loop* 

A liquid organic hydrocarbon will fulfill the requirements listed 
above j, but also will present some unsolved problems* It will (l) tend 

V 

to increase the size of the reactor by decreasing the hydrogen modera^ 
tion relative to a water cooled and moderated reactorg (2) tend to 
increase the reactor size due to its less efficient heat transfer 
capabilities compared with water and (3) exhibit both thermal and 
radiolyti© breakdown in which the original hydrocarbon willg primarily 
by polymerization^ change into longer chain hydrocarbons with the 
evolution of relatively large amounts of gas and a great increase in 
viscosity of the coolant^ 

The latter problems will be investigated at the National Reactor 
Testing Station on the Qrganic^Moderated Reactor Experiment (ORME)© 

This experiment is not a prototype for an organic power reactor but 
is being built to test the feasibility of its use as a reactor coolant* 
At the West Milton Annex of the Knolls Atomic Power Laboratory the 
Atonde Energy Commission plans to investigate the feasibilily of an 
organic moderated and cooled reactor for naval propulsion* This will 
be conducted on the Naval Organic Reactor Experimentg (NOEIE)© 
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The heat transfer properties of the various liquid organic 

hydrocarbons that are presently under consideration are not generally 

available* Some presently under investigation are diphenyl, 

0-terphenyl, M-terphenyl and P-terphenyl* A proprietary heat transfer 

medium ”Dowtherm" is available commercially and is composed of 

diphenyl oxide and diphenyl* It has a working range of 

7^0“F* at pressTires less than lh$ psig* Since ’’Doirtherm" closely 

resembles the lydrocarbons that will be investigated on the ORME and 

[ 2 ] 

its heat transfer properties are readily available it will be used 
in this investigation* 

The fuel element chosen for this investigation is of the MTR 

[3] 

(Materials Testing Reactor) type* An effective length of element 
was set at 1^8"* '‘‘his length was chosen so that the overall height 
of the reactor, control devices and shield would fit into the present 
day engineering spaces of a cruiser* Any additional height would have 
an adverse effect on stability, since it would result in a higher 
shield deck and thus an increase in the height of the center of 
gravity of the ship* 
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Fuel Element Dimensions 



t « 0.020" 
m 

t “ 0.015" 
c 

W » 2.5" 

L - 1|8.0" 

(Case 1) 
(Case 2) 
(Case 3) 




- 0.115" 

s 

= 0.0919" = 0.8 t^ 
s s 

3 1 

t « 0^069” = 0*6 t 

s s 



0 . 0 /^" 



PROCEDURE 



Heat Transfer ° 

The controlling variables in this study were (l) the maxiinum 
surface temperature of the fuel element or the maximum temperature 
at the centerline of the fuel element and (2) the maximum coolant 
speed. The maximum surface temperature of the fuel element may be 
controlling^ since the preliminary investigations on the ORME have 

[U ] 

shown that above 800® F the thermal breakdown is increased greatly. 
However, liquid loading reqtiirements of a naval vessel make storage 
of used coolant a simple matter and in the interest of higher 
efficiency some breakdown will be accepted. Therefore, a fuel element 
surface temperature less than 8^0® F will be accepted. If the maximum 
fuel temperature at the centerline is allowed to exceed a value of about 
1130® F, a phase transformation of uranium is encountered in which an 
anisotropic expansion and possible rupture of the fuel element may 
result. Therefore, a second design criterion will be required that 
the fuel element centerline temperatxire be loss than 1130® F. 

Preliminary studies on the ORME have also shown that coolant 
speeds greater than 15 feet per second result in excessive erosion. 

As a third design criterion, it will be required that coolant speed 
be less than 15 feet per second. 

The first investigation was with a reactor coolant outlet teit 5 >er- 
atvire of 750 ® F with various temperature rises across the reactor. 

This investigation was amplified by varying coolant speed while holding 
the temperature rise across the reactor constant. Each of these 
situations were further studied by considering the MTR fuel element 
modified to the extent that the plate spacing was narrowed to 80^ 
and 60^. 

• Nomenclature for Heat Transfer on pages 22 and 23. 
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Though it was known in advance that the speed would be kept at or 
less than 1^ feet per second calculations were pursued with higher 
values* This was done to show that a reduction in size can be accom- 
plished in the futurej if the erosion problem at high coolant speeds 
can be overcome* 

g 

A heat rate requirement of 6*06 x 10 Btu^r was chosen based on 
the full power steaming rate of a tjqsical cruiser* 

g 

With this value of Q (6«06 x 10 ) a coolant outlet temperature^ 
coolant temperature rise^ coolant speed and plate spacing are chosen* 
Since the inlet and outlet temperatures are known the average values 
of density and specific heat may be calculated* Formulae 1^ 2^ 3j> and 
h are used respectively to calculate w^ G and n* 



Q 

G AT 
P 



G A 



f 



w 

A » ■ ' 
f 3600pv 

w 

G * To » 

P 360QA 

f 

^f 

n ® ““ 

*,■ 

Since the dimensions of the ’:^el element for each case are known, 
both the hydraulic diameter and the total flow area can b!fr calculated* 
With the product of the mass flow rate and the hydraulic diameter the 
produce of the heat transfer coefficient and the hydraulic diameter can 
be read from the graph in the "Dowtherm” booklet* Formulae 6, 7s 
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and 8 are then used to calculate, respectively. Q , AT .«> 

^ Y fuel clad 

and ©• 



V V 



2 

t Q 

. m V 

AT « T - T = — *• 

fuel o 1 i|k^ 



t Q t 
m V c 

AT - T - T * 

clad o 2 2k 



Q t 

- _ „ V m 

© ■ T - T ■ — — ’ 

2 b 2h 



( 5 ) 

(6) 

( 7 ) 

( 8 ) 



Datum ftom the calculation of a coolant outlet temperature of 750® F 
for various values of temperatiure rise and speed for cases 1^ 2 and 3 
are in Tables Ip Up and HI and Figures Ip II, and HI respectively® 

The cladding temperature existing for each of these conditions is 
plotted on Figure IV, 

The narrowest plate spacing, case 3p gives a significantly smaller 
reactor® Since the pumping power of any of these cases proved to be 
very small compared xd.th the total power output, it was decided to 
investigate the result of a decrease of coolant outlet tenq)erature for 
this case only. Coolant outlet ten^ratiire was reduced to 700® V and the 
speed was held at 1^ feet per second. Datum from this calculation are 
in Table IV, These data were plotted as a dotted line on Figure IV • 

For these twelve conditions the fraction of uranium in the uranium- 
zirconium dispersed fuel plate was calculated and plotted on Figure V, 

It is desirable that the fuel be of such a composition that it contributes 
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TABLE I 



Datum for coolant outlet of 750® F 
and wide plate spacing 
(case 1) 



REACTOR AT * 200 ®F 








V 


G 2 


h 2 


Af 


n 


( ft/sec ) 


(Ib/sec-ft ) 


(Btu/hr-ft - 


®F) (ft2) (no. 


of plates) 


5 


236„1 


686 . 


5.36 


2^685 


15 


709.0 


I 6 II 0 


1.79 


896 


25 


1180.0 


23iiO. 


1.07 


538 


50 


2360.0 


ii090. 


o. 5 it 


269 


V 


«v == W®, 


T„-Ti T^-T2 


T2»T^ T 2 


^u 


(ft/sec) 


(Btu/iir=ft ) 


(of) (®f) 


(OF) (op) 


■ % 


15 


I,.86 


20.7 62.2 


215.8 1015 . 8 * 


25.lt 


REACTOR AT * 150®F 








V 


G 


h 


A- 


n 


(ft/sec) 


(ib/sec-ft^) 


(Btu/W“ft^- 


*F) (ft^) (no. 


of plates) 


5 


231.9 


682. 


7.26 


3^639 


15 


705.0 


1591. 


2„li2 


I 921 I 


25 


1160.0 


2320 . 


l.lt 5 


727 


50 


2318.0 


UOiiO. 


0.73 


361t 




Qv X 10® 


T.-T]^ Tj^-Tg 


T 2 “Tb h 


fu 


(ft/sec) 


(Btu/4a?-fV ) 


(*F) (®F) 


(op) (.p) 


% 


15 


3.60 


i5.it lt 6.0 


197.3 9lt7.3^ 


19.5 



excessive tenqperature 
'^rani\3m content too low 
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Table I (Continued) 



C, REACTOR AT » 100»F 








V 


G 


h 


Af 


n 


(ft/sec) 


2 

(lb/seC“ft ) 


2 

(Btu/W-ft - 


•F) (ft^) (no. 


. of plates) 




227o 


669 • 


11.01 


5i,520 


15 


681* 


1569. 


3.67 


1,839 


25 


1135. 


2295. 


2.02 


1,103 


50 


2270 


iiOOO. 


1.10 


552 


V 


JC lo8 


To-Ti T^-Tj 


^2”^b ’’2 


fu 


(ft/sec) 


(BtuAr-ft ) 


(OF) (*F) 


(op) (.p) 


% 


15 


2.37 


10.1 30.3 


133.8 883.8* 


13.5 


D, REACTOR AT * 50 ®F 








V 

(ft/se©) 


G 

(lb/sec“ft ) 


(Btu/hr-ft - 


®F) (ft^) (no. 


n 

of plates) 


5 


221.8 


659. 


22oli0 


11^230 


15 


667.0 


1515. 


77J*5 


3,730 


25 


1110.0 


2275. 


h»h7 


2^2UO 


50 


2218.0 


3910. 


2.2li 


1,123 


V 


«v == 


T.-T T -T. 


’’2“^ T 2 


f 


(ft/sec) 


(Btu/hr=ft ) 


(.F)-*- (op) 


(op) (op) 


U. 

% 


15 


1.17 


h.97 15.0 


66.2 816.2 


6.99 



* 

excessive temperature 

Y 

uranium content too low 
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TABLE n 



Datum for coolant outlet of 
coolant speed of 15 ft/sec, and medium plate spacing 

(case 2) 



-AT 


G 


h 2 




n 


X 10^ 


(*F) 


(lb/sec“ft^) 


(Btu/hr-ft ■ 


-•F) (no. 


of plates) 


(Btu/W-fu ) 


200. 


709. 


1671. 




1121. 


3.89 


150. 


705. 


1660. 




15L5. 


2.88 


100. 


681. 


1639. 




2300, 


1.89 


50. 


667. 


1631. 




1{660. 


0.9U 


AT 




VT2 


T -T 
,2 b 


T 

2 


^u 


(op) 


(op) 


(op)- 


(ep) 


(•F) 


% 


200 


16.6 


li9.7 


193.9 


9 ie. 9 l 


19.6 


150 


12.3 


36.9 


lhho2 


89li.O 


16.1 


100 


8.1 


2U.3 


96,3 


81i6.3 


11.0 


50 


lt.O 


12.0 


I47.6 


797.6 


5.6^ 



temperature excessive 
Y 

uranium content too low 



n 



TABLE ni 



Datum for coolant outlet of 7^0®?, 
coolant speed of 15 ft/sec^ and narrow plate spacing 

(case 3) 



AT 

(*F) 


G 

(ib/sec-ft ) 


(Btu/hr-ft - 


.®F) (no. 


n 

of plates) 


10®, 

(Btu/hr-ft^) 


200 


709 


1796 




1U96 


2.9 


150 


705 


1775 




2020 


2.2 


100 


681 


1731 




3061 


l.U 


50 


667 


1718 




6230 


0.70 


AT 


T =T 


T ~T 
1 72 




'^2 


f 

u 


(®F) 


(»F) 


(«F) 


(«F) 


(«F) 


% 


200 


12,li. 


37.3 


135.2 


885.2* 


16.3 


150 


9.2 


27.6 


101.3 


851.3 


12oli 


100 


, 6.1 


18.2 


68.5 


818.5 


8.1;^ 


50 


3.0 


9.0 


33.9 


783.9 


l»o3^ 



temperature excessive 
Y 

uranium content too low 
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‘ Coolant speeid 3a ft/S't© 
13 



nraibtG* of platfs 






1 




it/ soe 



niaaibe? of plat*s 



liMi ••••■•■»•(• iiaaa ••• I 






;ss:i2:xssu»: 



txmxxmxizxR 

MMM 



Ml 

ixxxxxxixx: 








mmm 




TABLE IV. . 



Datura for coolant outlet of 700*F, 
coolant speed of 15 ft/seCj and narrow plate spacing 

(case 3) 





G 


h 


(*F) 


(ib/sec-ft^) 


(Btu/hr-ft 


200 


739.5 


1826 


150 


72608 


1812 


100 


712.5 


1775 


AT 


T -T 
® 1 


T -T 
1 2 


(®F) 


(op) 


(•F) 


200 


21o3 


6U.0 


150 


15.9 


hi. 9 


100 


10.5 


31.5 





n 


X 10® 


’F) (noo 


of plates) 


(Btu/hr-ft^) 




871 


5.0 




1164 


3.7 




1767 


2.5 




'^2 


f 

It 


(•F) 


(•F) 


LI 

% 


228.U 


924 . 4 * 


26.0 


172.2 


872.2 


20o3 


115.8 


815.8 


14.0 



teraperature excessive 
Y 

uranium content too low 
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FIGDEE V 








to th* structural strength of the reactor. A uraniinn content of 10^ by 
weight or greater has desirable strength characteristics. 

Of the six conditions that have a fuel element surface temperature 
of less than 850® F only two have a uranixm fraction of above lOjC. 

These are shown in Table V, 



t 

CO 


T 


case 


n 


^2 


750 


100 


3 


3061 


818 


700 


100 


3 


1767 


815 



Table V Particularly Interesting Cases 

Figure VI shows the variation in number of plates for case 3 
holding the coolant speed at 15 feet per second and varying the 
temperature rise across the reactor. 

The size resulting from the condition with a coolant outlet 
temperature of 700® F is better for the following reasons# 

1) it is by far the smaller, 

2) the increase of uranium fraction by 2 1/2^ gives a 
significant increase in the strength of the fuel region. 

With this condition the fuel element surface temperature will be 
less than the prescribed value of 850® F. The amount of damage that will 
result to the coolant from this temperatxire and the accompanying 
radiation is uncertain. However, for shipboard applications it does 
not appear to be a serious penalty. 

In the present day cruiser design for torpedo protection there is 
available a large amount of bunker fuel oil to partially fill the 
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lasners of the torpedo protective system* Certain layers of the torpedo 
protective system are required to be liquid loaded* This system could 
be used for the storage of the organic moderator coolant* The unradiated 
coolant could be stored in the inner layers of the system and that 
irtiich had been irradiated would be stored in the outer layers to protect 
the crew from any induced radioactivity* The long chain hydrocarbons 
wDxild tend to settle to the lower parts of the tanks where they co\ild 
be removed during periods of overhaul* With the installation of a 
purification system this process of separation might be handled con- 
tinuously aboard ship in the same manner that lubrication oil is 
purified and recycled* 
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Nomenclatiire 



f 

n 



P 

D» 

f 

u 

G 

h 



flow total area, (ft ) 

flow area per plate case n, (ft^plate) 

heat transfer area, (ft^plate) 

average specific heat, (Btu/lb ®F) 

hydraulic diameter, (inches) 

fraction of uraniinn in fuel region 

2 

mass flux, (ib/sec-ft ) 

2 

film coefficient, (Btu/iir-ft -®F) 
thermal conductivity of cladding, (Btu/hr-ft-*F) 
thermal conductivity of cladding, (Btu/)ir-ft-®F) 
fuel length, (inches) 



n 



number of fuel plates 



CO 



m 



total heat generation rate, (Btu/fu ) 

3 

heat generation rate per unit volume of fuel, (Btu/hr-ft ) 

fuel cladding thickness, (inches) 

coolant outlet temper at\u*e, (®F) 

fuel thickness, (inches) 

spacing between plates, (inches) 

bulk temperature of coolant, (®F) 
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Nomenclature (Continued) 



fuel centerline temperature, (®F) 
fuel-clad interface temperature, (®F) 
fuel element surface temperature, (®F) 

2 

U overall heat transfer coefficient, (Btu/hr-ft -®F) 

V coolant speed, (ft/sec) 

w mass rate of flow, (ib/hr) 

W fuel width, (inches) 

AT overall temperatTire difference across reactor, (*F) 

^ average density, (ib/ft ) 

6 film temperature drop, (®F) 
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Batic R»actor Phyaics 



Criticality calculations were made using a one group reflected 
homogenous model with case 3 9 narrow fuel elements* Based on heat 
transfer considerations, inlet and outlet coolant ten5)erature were 
set at 600* and 700* F, respectively, with a coolant speed of 1^ ft/sec* 
The assun^tions in this calcvilatlon were 

1* The structural members inside the reactor were of 
type 3h7 s.tainless steel and composed ten percent of 
the total reactor wei^t, 

2* A fast fission factor of unity, 

^ 235 

3* The uranium of fuel is enriched to 90% in U , 
li* A right circular cylindrical flux distribution of a cosine 
axially and a Bessel function radially, 

5* A beryllium reflector of 236 cm placed around the reactor* 
The reactor was treated as a homogeneous medium* By a one group 
calculation the volumes of each constituent of the reactor were cal- 
culated from the dimensions given previously in the heat transfer sec- 
tion* These volumes were then converted to weights and the volume 
fraction of the total reactor as shown in Table I* From the cross 
sections given in BML-32^^ ^ ^ effective cross sections were calculated 
for the operating temperature of the material* With these cross 
sections the reactor parameters as shown in Table II were calculated* 
These give an effective multiplication factor of 0*89li* This value is 
sufficiently close to tmity that final ad;justinents are made by the 
computing machine* 

The primary purpose of the calculation is to obtain an approximate 



valu* of th« critical radius which is needed as an input to the 
computer solution* 



TahXe I 

Reactor Materials 



Volume 
(cm ) 

Fuel 6*95 X 10^ 

U 
Zr 

Clad ( Zr) 15.01 x 10^ 

k 

Stainless Steel 1.99 x 10 

Cr ly weight) 

Ni (13^ ty weight) 

Fe (Tl^ty weight) 



Mass 

(Kgm) 


Volume fraction 


li 95 . - 




69.5 


q'*oo68 


li 25*5 


^ 0*li08 


966 




160 




28 *li 


0.007 


17.35 


o.qoU 


112 . 


0.027 



Coolaijit 26.3 x 10^ 219*5 

[73*5 Diphenyl oxide by yeight (C H OC H ) and 

6 5 6 5 

26.5^ Diphenyl by weight (C,H ) ] 

652 



C (06*91 by weight) 
H (06.03 by wei^t) 

0 ( 86*59 by weight) 
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0*U70li3 

0*03292 

0*377lt 



Tabl« II 



Reactor Parametera 



tr 



Berylliun shield thickness 

5 ^ ) 

- 1 



0 . 






1*3 

tr 1-3 



33.7D 



2 V 
». f 

X 



2 O 
L : 



p(E) - «p 1 10-^'* 

' ”28 






■V 



^a( total) 



= /^fp 6 

g M i 



r 



E 

“it 



B 



2,U0g 

R + d 



H+26z 



3 = albedo of Beryili\im 

Ir iSssS— — 

ef f “ 1 +b2l2 



23*6 cm. 

0. 28331t cm"^ 
1.176!tU era 

39.599 cm 

0;i208 cra"^ 

9,7387 cm^ 

0.726 

1.9692 

1. U296 
o;17598 

2 

3li3 cm 

1.332 X 10“^ 

0.89U 



2 

cm 
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Application of the Flux Equations 



The S approximation is a method of differencing lii which the 
n 

independent variable is quantized into n intervals, and intermediate 

'values of a function of the -variable are approxima-bed by a straight 

line be-tween points at the boundaries of the inteiVal* This method 

[61 

of numerical integration was developed by Bengt G* Carlson • 

The -vector flux at a point within a reactor can bo ob-bained 
the use of the me-fchod* The equations necessary are deri-ved in 
Appendix -A. for -two cases of neutron velocity? 

Case A - The component of the vector flux in the axial direction 
is grea-ter than i-ts component in -the radial direction* 

Case B - The component of the vector flux in -the radial direction 
is greater than that in the axial direction* 

These two cases are illus-ferated in Figure I using a grid which is 
used to difference -the -transport equation in the r and z coordina-bes, 
r being the radial coordinate and z the axial coordina-te* 

Figure I 

The -two cases of flux at a grid point 
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Th« general flux equation for case A ist 




where h Is a value which is one greater or less than j, depending 
upon the direction of the radial conponent of vector flux, and n is 
a value which is one greater or one less than n, depending upon the 
direction of the axial con^jonent of vector flux* Definitions of the 
other quantities can be found in the nomenclature part of this 
section* 

For (Jase B the flux at a point is defined tgr the following 




where A for Equation (l) and (2) is defined bya 



A 















k 




8,n 



V'" 

b,m 





( 3 ) 
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It should be pointed out that the total macroscopic cross 

section^ ^ is evaluated using the material properties at the 

center of the position grid of Figure I, but are designated by the 

grid values at the corner diagonally opposite the point 

i»e.| the point h^n« Further^ r is actually the mean radius and is 

h 

r * r 

equal to 3 h • 

2 



The source at a point is defined by the following equations 

I V * Hvy., 

k*l f 3 ,m k»l ^ 
where the scalar flux 0 is equal to 



i H I P 
f -0 f i-o i 









The fission cross section and transfer probabilities^TT^ are 
evaltiated using the material properties at the grid point at which 
the flux is being calculated* The data necessary for the solution 
of the flux at a point is summarized in Tables X and H/. 



Table 11 



Data^f or the approximation 



i 


J. 


JL 


i 


i 


i 


i 


0 


-1 


-1 


0 




1 


•K).3183 


1 


0 


-1/3 


V3 


-2/3 


2 


•fO,361(3 


2 


1 


2/3 


k/3 


1/3 


2 


■K)*3183 
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Table IE. 



f 

0 

1 

2 

3 



Data for Gauss quadrative approxination 



-0,8611363 

-0,33?9801 

+0.3399801 

+0,8611363 



-1.69Uli625 

-0,3197281* 

♦0,3197281* 

+1.691*1*625 



+1.9677051* 

♦1^06331*0 

+1.06331*0 

♦1.9677051* 



+0,91*7851^ 

♦0.65211*52 

♦0.65211*52 

+0,31*7851*8 



Not all pairs of values for i and f are rlilevant for each of the 
cases A and 6 as a consequence of the defining inequalities for these 
cases. The pairs of relevant values are discussed in connection >d.th 
Figure II. A general point may be defined to be any grid point of the 
reactor other than a point at the center, on the centerline or in the 
mi^plane of the reactor. A point at these other places is called an 
oidgin point, a center^JLin* point or a midplane point, respectively. 
For Case A for ^rtiich|b.| <|h^|and in the case of the first quadrant 
flux, there is a further requirement (Table I, Appendix A) that b. 
and have negative signs. In this quadrant there are three Case A 
f^Ltoc equations, one for each of the pairs (i*0, f-0), (i*l,f«0), and 
(i«l,f*l). Similarly, for Case B, the requirementjbj >|h^|lead|| to 
only one relevant combination. 

The applicable combinations and those for the remaining three 
quadrants are summarized in Table HI, 
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Figore II 

Flux Configuration 
for a particTilar energy group 
at a general grid point 
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Table IH 



Allowable combinations of i and f 



Quadrant Case A Case B 



1 


i - 0, f - 0 
i - 1, f - 0 
i - 1, f - 1 


i - 0, f = 1 




i « 2, f - 0 


i - 2, f - 1 


^rd 


i = 2, f » 3 


i - 2, f *2 




i - 0, f - 3 
i = 1, f - 3 
i = 1, f - 2 


i « 0, f • 2 

f 



The configuration of the flux at a grid point which results from the 
allowable combinations of i and f are shown in Figure II* If the flux 
is assumed to be distributed in^+1 energy groups there will be 12 (J^ 
flux equation associated with every grid point* For the case considered^ 
there are three energy groups and, therefore, thirty-six flux equations 
for each point* 

[7] 

The following boundary conditions are assumed* 

1* Continuity of i^ux exists at the centerline and at the 

midplane oX the reactor, and further the flvx is symmetrical 
about these boundaries* This permits the calculation of 
the flux in only one quarter of the reactor* Figure III 
shows the flux distribution at these various special points* 
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Figure HI 




2« A second boundai*y condition is that at the outside 
boundaries of the rea,ctor there exists no Inward flux* 

In order to satisfy this requirements, an extra set of grid points is 
placed just outside these boundaxdes, and at these points the vector 
flux whose argument describes inwardly moving neutrons is asstmed to be 
zero* Note, however, that leakage does exist at these fictitious 
points and at the outside boundary points of the reactor* 
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Nomenclature 





P 

i 







Values obtained by dlfferencjbgnlnd'lntegratloa^ of 
the transpoi*t equation overO*J f and defined by 
Equation (?4 to 7d) of Appendix A 

Valume depending upon the quantized values of |x and 
defined by Equation (?e and 7f)of Appendix A, 

Values obtained by the use of the Gauss quadrature 
approximation to sxm the flux over (i^paaqi^e. 
defined by Equation (13) Appendix A» 

Values obtained by the use of the approximation 
to sun the flux over ^ space, defined by Equation (12) 
Appendix A* 

mean radius between two adjacent grid points 
Neutron source 

Fraction of fission neutrons which are produced in a 
particular energy group 

Number of neutrons produced per fission* 



3k 




The probability per unit length of neutron travel 
that a neutron is scattered in group k and that a neutt^ 
emerges in group A, from the scattezang event* 



Macroscopic fission cross section 
Total Macroscopic cross section 



neutron flux i^ch is classified according to speedy 
Tf colatitude and azimuthal direction^ cosines p and »» 
and position r i(nd z* 



A 



A 

m 







UJ 

uu' 



spacing of grid in the radial direction* 

spacing of grid in the axial direction* 

cosine of the an^e 0 between the axis and the 
direction in which neutrons move* See Figure I 
Appendix A* 

Cosine of the angle y/ * See Figure I Appendix A* 
See defining equation^ page A-12 Appendix A* 

See defining equation, page A -17 Appendix A* 
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Nomenclaturo 

Superscripts and srabscripts 



i index labeling the quantized values of , 

f index labeling the quantized values of |ju 

j grid position index in the radial direction. 

\ 

Bi grid position index in the axial direction. 

^ neutron speed index. 

k neutron speed index of group in ^ihich neutron is 

scattered • 

n axial grid position diagonally opposite grid point 

at which flux is being calculated. 

h radial grid position. 



36 



Coaiputer Introduction 

A computer program is a sequence of logical instructions telling 
the computer what mathematical operations are to be performed on the 
data. Therefore, there are two groups of inputs, the data and the 
instructions. The data fed into the computer are stored in the fast, 
core memory if there is sufficient room. When the space is not 
adequate, some of the data are stored in magnetic drums. T^e only 
diff^ence in the two storage processes is that the data in the 
drum takes longer to locate and bring into the core memory# Each 
number or piece of data is stored in a register or cell the location of 
which is designated by a number called the address of that register. 

The instructions are stored in sequence and are performed in the order 
in idiich ttiey are stored unless the program directs otherwise# No 
unused registers are permitted to interrupt the storage as these will 
stop the computer# 

The instructions given to the computer consist of two parts, the 
operation and the address sections# The operation section is first 
and uses letters to tell the computer the mathematical operation to 
be performed# The address section gives the address of the register 
containing the data on which the operation is to be performed# The 
instnictions recognized by Whirlwind I are listed in the glossary# 

For Whirlwind I there are many useful devices that make the 
programming ranch simpler# ^ of these devices used in this 
program are floating addresses, preset parameters, and subroutines# 

Probably the most useful device is the floating address. This 
floating address or "flad” is used to eliminate the necessity of 
referring to the absolute register address in the core memory when 
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Tndting the program* The absolute register address of any part of 
the program must be known primarily so that the initial or final 
instruction of data is not placed in a register containing part of 
the subroutines fixed in the memory of the computer* 

the program key instructions on which mathematical operations 
are to be performed or which will be referred to later in the program 
are given floating addresses* These addresses consist of a letter and 
a number^ for example al, a2, etc* letter except o and ell 

may be used and the sum of all the numbers following these letters may 
not exceed 255 ® When the program is read into the computer, it will 
assign absolute addresses or register numbers to these flads* A 
table of the floating addresses is produced by the computer with the 
solution unless it is specifically suppressed in the computer* ^his 
table gives the absolute address of every assigned floating address 
and is particularly helpful in ^•^ouble shooting a faulty program* 

Another useful device is the preset parameter* ^ese are 
numbers other than data that are used in the program* Preset 
parameters are frequently used to designate the number of times that 
a calculation is to be performed* Preset parameters must be assigned 
a value before they are used, and they will retain this value until it 
is specifically changed* Preset parameters are designated by two letters 
and a number* ^e first letter must be a p, u or z and the second 
letter and the number may be ai^hing, except the letters o and el may 
not be used* 

A third useful device is the subroutinp* Its usefulness comes 
from the fact that it is a conqjlete entity in itself and may be construc= 
ted and trouble shot independently* ^or the build up of a complete 
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program it is possible to start with one of the most basic mathe” 
matical operational sequences that will be performed repea tedly, 
program it, perform the hand calculations of the sequence and 
then place the routine on the computer* Modification of the 
subroutine or correction of the hand calculations must continue 
until the results of both are in agreement* 

In arranging these subroutines in the main program ex-bTeme care 
must be used* The main program must be interrupted to call the sub- 
routine into use and after the subroutine has performed its function 
it must direct the con 5 >uter back to the proper place in the main program* 
An exan^jle of a subroutine is as follows* 

Main program 

o 

o 

o 

isu p 10 

tl^ isp X 1; take the next instruction from register 

X It and continue from there 

imr a 5 forms product of (=m2) and (a5) 



i STOP 

xhf ita X 5 

ica m 2 1 
its a 1] 



end of program 

transfer into the address section of the 
instruction in register x5j the address 
that is one more than the address tl* 

places m 2 in a 7 
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ics m 2 
its a 8 

x5s isp 0 returns to main program one register 

after t Ij i»e»j, register containing 
imr a 

Since the numbers used in this program were over four decimal digits 
in lengthy floating point arithmetic was required* In this system two 
registers are required to express a number* Each of these registers is 
capable of containing 16 bits of binary number information* The first 
bit of the first register expresses the sign of the number and the 
first bit of the second register expresses the sign of the exponent of 
2* Six bits are required to locate the binary point of the number 
leaving 2U bits to express the actual number* In certain stages of the 
program it is desirable to handle single length register numbers* The 
instructions for the handling of these numbers are similar to those for 
double length registers except that the i does not precede the instruc- 
tion* When operating on double length registers the computer is said 
to be in the interpretive mode, when on single length registers it is in 
the Whirlwind mode* To enter the interpretive mode from the Whirlwind 
mode it is necessary to use the instruction "IN” and to reverse the 
direction it is necessary to use the instruction "OUT"* 



places -m 2 in a 8 



GLOSSARY OF INSTRUCTIONS 



TERM 



DEFINITION 



ica X 
ics X 
iad X 
isu X 
iinr X 

idv X 
its X 

iex X 
Isp X 



clear MRA and add contents of 
register x* 

clear MRA and subtract contents 
of register x* 

add contents of register x to 
idiat is in MRA* 

subtract contents of registOT x 
from what is in MRA* 

multiply what is in MRA by the 
contents of register x and round 
off product to fifteen digits* 

di'vide contents of MRA by the 
contents of register x« 

transfer the contents of the MRA 
to register x losing what was 
previously in register x* The 
MRA remains unchanged* 

exchange the contents of the MRA 
with the contents of register x* 

take the next instruction from 
register x* Does not affect the 
contents of the MRA* An uncon= 
ditional instruction used to break 
the sequence of operations* 
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Glossary of Instructions (Continued) 



TERM DEFINITION 

icp X if the contents of the MRA are 

negative^ take the next instruc- 
tion from register x and continue 
from thereo If the contents of 
the MRA, are positive, take the 
next instruction in sequence* A 
conditional instruction used to 
break the sequence of operations* 

isc X select co-unter number x* Without 

this instmction counter zero 
will be used whenever +C appears* 

icr X cycle reset® Sets index register 

of counter to zero and criterion 
register to x® 

ict X cycle transfer* Increases con- 

tents of index register by one® 

If contents of index register greater 
or equal to contents of criterion 
register, set index register to 
zero and do next instruction 
in sequence* If contents of 
index register less than contents 
of criterion register, take next 
instruction from register x® 

ici X increase the contents of index 

register by number x® 

icd X decrease the contents of index 

register by number x. 
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Glossary of Instructions (Continaed) 



TERM 



DEFINITION 



iti X 



iat X 



ita X 



iTOA 

iMQA. 

iSOA 

IFOR 



transfer the right eleven digits 
of the index register into the 
right eleven digits of register x« 

add contents of index register to 
the contents of register x and 
store the result in the index 
register and register x* 



replace the address section of the 
instruction in register x with the 
address that is one more than the 
address of the register containing 
the last isp or icp if the contents 
of the MRA are negative© 



record the contents of MRA on 
direct printer (typewriter)© 

record the contents of MRA on 
magnetic tape for delayed printing* 

record the contents of MRA on 
oscilloscope for photographing* 

this instiruction provides an 
automatic device for obtaining a 
suitable layout of output data 
in colximns, lines or blocks* 



NoteJ The above output instructions are usually followed by a series 
of letters and numbers that indicate the desired form and 
arrangement of the output® 
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Compater Logic 



The jjrogrem starts at 1 in the main program by entering the 
data, program and preset parameters. The initial calculations are 
with the reactor radius equal to the value placed in register rl*. 

With this radius the program then shifts to the flux iteration and 
convergence test, ii. In this routine parameters are set to control 
the number of iterations allowed at one radius and to traverse the 
conpiter thru all points in the reactor. With a value of j and m, 
defining a point in the reactor, the program shifts to Flux, U.5> 
where the addresses on the drum of the required data are calc\ilated, 
the data entered in the .fast memory and the routine appropriate to 
the point is selected and control is shifted to that routine. 

The first calculation at any point is the scalar fli^ for all 
throe groups and the second is the input from fission and scattering 
from an upper energy group. The first point to be calculated in each 
cycle is the origin. This is the point of normalization of the 
inputs to the fast group. The normalization constantj,^, is defined 
as the reciprocal of the fast group input at tdie origin. With this 
constant all other inputs (all groups and all points) are normalized. 

As explained in the previous section of the procedure certain 
quadrants of the flux are required for each point depending on its 
location. Prior to sequencing to another point all fluxes appropriate 
are calculated. 
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fast sotirce ii^ut is then eoBQ>ared to the preTlotw fast 
source input* If the new value is not within a set specified fractional 
value of the previous value a negative quantil^ is Introduced into s 10* 
When the entire grid has traversed a check is made to see if the 
maxunum allowable attempts at source convergence have been made* If this 
maximum number has been made the counter will stop* If not a check 
will be made of switch s 1(^ As indicated in routine kt another eyele 
will be made thru the reactor if « 10 i» negative^and if ^ 10 is positive 
control will return to the main program to repeat the above procedure 
based on the radius placed im r 3« 

f 

lifhen a gamma is calculated for this radius^ control is again 
returned to the main program but at a different point* All subsequent 
returns are at this point* A check is made to determine if the 
maxlm^ allowable nomiber of Iterations of radii have been attei^tedo 
A test for critlcali'ty is made to see if gamma is within a 
specified tolerance of unity* If it is not^ a new radius is found by 
extrapolation or interpolation based on the most recent values of R and 
y and the best previous value of R and » In this determination of 
radius^ the gamm a that is closer to tuiily and the associated radius arm 
preswved for use in later extrapolation* The previously outlined cal- 
culations are repeated with this radius* If ^ lies within the specified 
tolerance of unity control remains in the main routine and the print-opt 
i^tine^ which has .'tym stored in the drum^ is entered into the fast 



memory* 



Th« computer is then sequenced thru the grid calculating and 
printing out scalar fluxes at each point for each group. 

If the flux pattern is not as desired the properties of the regions 
are changed and a new tape of data cut. The program is then run again to 
determine the flux pattern. 







U7 





on subsequent cycles 



1»8 







h.$ 



Ebd-t to 

Grid traverse (lio6) 



Ggid Point Sensing 
Enter flrom Address calculation and 
data entey (it^ol) or calculate drum 
address of required data and enter 
daU (8.1) 



lSa2 






\S2^1 



c 






>^--c 



< 



\ 62»2 



C 



mSO 












mSO 



Ebd.t to 



> 



> 



< 



y 

Origin 

Program 

Ot5.3) 

or 

Ori^n 

Scalar 

Flux 

(i»53.l) 



< 






V 


Y 


1 


Centerline 


Midplane 


Interior 


Progjcam 


Program 


Program 


(li5.5) 


(li5.li) 


(Ii5o6) 


or 


or 


or 


Centerline 


Midplane 


Interior 


Scalar 


Scalar 


Scalar 


Flux 


Flux 


Flux 


(i»55a) 


(ii5iiol) 


(156.1) 
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Qrigiln program 
Enter from Grid Point 
Sensing (lt^o2)5 origin 
exit 



li5o3 



U53ol 



U53o2 



l*53o3 



h$3ck 



k^3o^ 




ft“int (8.2) if 

|y-i s£f 



Exit to 

Throw source con= 
mergence switch 

(ii5o7) 
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Midplana program 
Enter from Gbrdd Point 
Sensing (hSo2)^ midplane 
exit 






l*53o2 



k$3oh 



kB3o^ 







Print (8*2) if 



Exit to 

Throw source con<= 
vergence switch 
(It5.7) 
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Centerline program 
Enter fr> 0 m Grid Point 

Sensing (li5o2)^ centerline 

e3d.t 




Print (8*2) if 

IM * fy 



Erlt to 

Throw sourbe eon^ 
vergence swits^i 

(hB.7) 



^3 



U56.1 



Interior point program 

rMiiMi ^ (^Trr-Wf q'L fci COMOMdCMaMaD 

Entar from Grid Point 
sensing ihBo2)g general 
point exit 






h^3oh 






h^hol 



h ^.2 



hSS.2 




switch {hS»7^ 



$h 





Tta>ow sougca convergence switch 
Enter fipom origins inidplane, 
centerline or interior point 
programs (li5f,®3 to ii5*6) 




Place flux in drum 



Grdd Tra vet-88' lto6 

Eater from Flux (fi«5) 
or rWnt (8.2) 




Exit to 

^ iteration check (li«7) 
or Stop 



Exit to 
Flux (lio5) 
Drum address 
calct^tion 

(8;d 



R?dnt out routine 
Enter from place pirint 



8 



out routine in MRA (7) 




Stop 
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RESULTS 



Niimerous runs of subroutines and of the entire program were made 
on Whirlwind with the program contained in Appendix This program 
was adapted for a reactor radixxs of 75 suid 95 centimeters and a height 

of 168 centimeters# The grid placed on the reactor was a 6 by 6 grid 

gi-ving for the first trial radius a grid size '15 cm by cm* 

The final run had programmed a first trial radius of 75 cm}' the 

maximum allowable number of iterations wi'th a cons'bant radius was 
•twenty and the tolerance of convergence of the fast source input was 
•three percent* The, computer stopped, as programmed, i: 5 )on conqjletion 
of "twenty iterations at a set radi'us ■when the fast neu-tron input at . 
all points "was not -wi-th the tolerance specified* The con^juter running 
time per i"teration is 'two minutes* The compu"ted values of scalar 
flux at each point and energy g|?oup w«re printed out for each 
iteration* The normalizing cons"tant decreased from 0*27 for the 
first iteration "to 0*01 for the twentie-th* Some "values of "the scalar 
flux at the following point were negative 

f 

a) interface of core and reflector, 

b) reflector, 

c) outside boundary of reflector* 
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DISCUSSION OF RESULTS 



The negative values computed for scalar fluxes are known to be 
unrealistic* The possible causes of this error are 

1) iJ 5 )roper programming of the theoretically derived 
flux equation, 

^ 2) grid size being excessive* 

Tihe theoretical backgrotmd for "Uie flux equations is sound, and 
these equations have been programmed for other geometrical arrang«nents 

Cn]' 

on other computefe and give correct flu± distributions * The 
need for a conqjlete hand calculation is well realized, but due to lack 
of time these hand calculations were not completed fully* Therefore, 
there is no reason to assume or expect that these equations have been 
programmed correctly* To accon^lish the required hand calculations, 
conqjuter verification of the hand calculations, and program trouble 
shooting it is estimated that another three to four weeks fixLl time 
would be required* 

For accuracy in evaluation^ of rieutron transport equations, the 
grid spacing should be less than or approximately equal to the mean 
free path of the neutron* If this condition does not prevail, the 

k - 

accuracy and reliability of the calculations decrease rapidly* 

A two fold increase in the nuidber of grid points will increase the 
time of computation by a factor of eight* Any increase in the time 
required for one iteration thru the reactor is not allowable since ndir 
two minutes are required* 
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The code developed is in the ‘interpretative* or ‘floating 
point* mode for Vftiirlwind^ for the reason that this mode is simple 
to learn and to trouble shoot* \ftiirlwind has also available a 
•^^lirluind' or ‘fixed point* mode in which the speed is increased 
by a factor of rou^ily 35* Use of the fixed point mode would 
make the two dimensional 

2 

practical* 



calculation on the Whirlwind coB^mter 



c<ajca:«sioNS 



1) The code as presently prograimned is incorreet* 

2) The grid spacing as used is much too coarse* This 
difii^culty is not inherent in the coded program* 

3) The mode of conputar operation is unsuitable for an 
iterative process such as this because it results in 
a eode that progresses too slowly* 
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RECOMMENDAnONS 



1) Prior to the use of computer a formal or semi^fonnAl 
coxirse in programming should be ccra 5 )leted* 

2) A grid spacing comparable to the mean free path of the 
neutron should be selected* 

3) The fastest possible mode of calculation on the computer 
available should be used* 

It) For a coded program of this difficulty a faster machine 
with a much larger fast memory is very desirable* 
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Appendix A 

Development of Equations Using Sg Approximation to the 
Boltzman Equation for Cylindrical Geometry 
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Equations will be developed for the time independent case using 



cylindrical geometry. 

The vector flux at a point in space is a function of two variables^ 
one describes distance along the vertical axis of the cylinder, and the 
other describes radial distances. ®ie axial variable id.ll be designated 
as p. and the radial variable asr^. Prom Figure I it can be seen thati 



» cos V 
p = cos 9 and 
sin0® 



{2 



Figure I 

Components of neutron flux 
at a space point. 



sin f » yl - r| 

A- 

and where -/X represents the unit 
vector in which direction neutrons 
travel at a point. Therefore, 

V 0 = i + i — - + i 



r i r 



and -/a ® i sin 9 cos 
r 

is 

i 



z ^z 



♦ A 




0 sin 0 sin cejs © 

ThereforAj, 

s sin 0 cos 



A 

-A. 



d r 



sin Q sin 
r 






cos 9 



11 
S 2 



If we substitute the expression for and p which were derived from 
Figure I in equation I 3 , we find 



A 

jn 



'a 






3 z 



( 1 ) 



dr r 

The quantity ® 0 represents the leakage of neutrons out of an 

element of volume per unit time in the direction per unit solid angle 



( 2 ) 



I 







In the J;ime independent case this leakage pltts the rate of loss of neutrons 
y. per unit time must equal the source strength^ So In problems 



in -which cross sections vary with energy and neutrons alter their speeds^ 

[ 6 ] 

the multi group approximation of the Boltzmann equation is introduced* 



or r 






•VV 



s 



(r^z^^ ) in Figure II can be 
represen-bed by the following 
equation where i represents 
the inter-val numbero 



7? Mesh 



(3) 



X y) 

where 2 is the total cross section within the energy group UC and 

a 

S is the neutron source term which results from fission and energy 

degradation of the total vector flux au the point* Slowing down is 

incorporated into the source term* 

Equation 3 is identical in form to the time dependent multigroup 

equation for spherical geometry in which the z term is replaced by a 

term representing the rate of change of flrix -with respect to time* 

The S approximation is a method of e-valuating integrals o-ver a 
n 

-variable in which the variable is quantized into n intervals and 
intermediate values of a function of the -variable are approximated by 
a straight line between points at the bo-undaries of the intervals* 

We shall consider ^ space as the first -variable which -will be into- 
grated* The flux at any point Fif 
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(ii) 



Fot the S^ approxiination i = +1^, 0^ and ->lo If we d«QP^ *^i” *|i 1 



as and ^CraZ^M.^*^ ^) as ^ we obtain 






i i 

If we substitute this in equation 3j we find 



~1 1-1* '5V /l<ai~1> ) H'li-i* // 

A^ ^ r A^ r r 

A >1 1-n^ A^ > z A 3 z 



Jl 

Q S 



■2( 



1.“ li.=l 






, £1 £ . , J.=l 

.. ^) V ^ ‘\i°" 1 



Since the sourcSg gj, is isotropic^ integrating this expression over Y) 

\ 

from |(|^ ^ ^ "fco we obtain n equations in the form 



£ 



Jtj. 



^ J.=l 



i9r ir l>z ^ ^ ± dr 

4 r-T Tr-r r- ^ g SjO = « gS« 

az ^ 



£ d=l 



i r 



(6) 



where 






(7a) 



'i^3A ^^“*\i“^i “'ll"! “ 
i '• 
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(7b) 



d. 






1-1 



3 



« * 2 
i 



(7c) 

(7d) 



U « ^ 

“\|1 - 



(7e) 

(7f) 



A different method of finite differencing is nsed for the other 
directional variables<^a since use of an expression analogous to (k) 
for (i in the Integration of the equation ($) ^th rjfspect to p would 
result in an extremely complex expression. Instead, the equation (5) 
is directly quantized with respect to p giving the resul1^»\ 

r i r f 0 z f 



Jj.-l«f 



( 8 ) 



\ “a"T” " i \ "#7~" ^ ^f ^ " *i 



SfS 



We may now proceed to eTaluate Uie sourfe integral at a point© 
For isotropic scattering and three energy groups of neutrons 

S(rj,z) « T^p ^2^ V(rpz) + V VCi^az)^ 



where 



k 



^(rgz) 



2tr 



♦1 +1 

i .{ 



i'X 



As in equation (li) J^(rj,Zj,p,^) is approximated by 

A ' T_1 

^(*’s2,p, ) ® ""“7- p (r,a,p) ♦ 



^ (raz,p). 



(9) 




If values of . ■ +lp 0 and =1 (S^ approximation) are substituted in 
the above expressior^we find that 

“ .1 P. s (12) 

d\t 1»U 1 

wher® P s 0 o 31835 P. ^ Oo363iia and P^ = 0„3l83 . 
o 1 2 

[5] 

By us® of the Gauss quadrattve method the remaining integral 
can be evaluated in the form of a sum of the products of coefficient 
times the vector flux where values are the n+2 positive roots of 



The table below lists the values of 



and the coefficient H used 

f 



f 

0 

1 

2 

3 






and the source term equals 





H 

f 


-008611363 


0®3ii7851»8 


“0® 3399810 


0®6521U52 


•H)® 3399810 


0„6?21ii52 


•K)® 8611363 


0®3lt785U8 


H = H the 
1 2 


scalar flux at j 







s - 1 '% ''s, (r,2) * ,J, ^ V . 



k=l 



(13) 



(Ih) 



where 

p is the propability that a fission neutron will have an energy^ and 
Ac 

T^s the probability that a neutron having an energy k will be 
scattered into the energy group <> 



The scattering probability requires further explanation® We define 

tji, 

as the probability that a neutron which is scattered in the k 

ith 



group will land in the H groupp k as the energy group in which tiie 

scattering occurs (l for fast^ 2 for epithermal and 3 for thermal) and 

A- as the energy group into which the neutron is scattered® THie 

energy range covered by scattered neutrons is E-oE^ and the proba- 

AE 

bility of landing in the interval AE is ® Let E be the energy 



E-aE 



of the neutron before scattering® 



E < E << E o 
k+1 k 






The neutron flux per unit energy interval is 



dE • 



Then dE 2 (E) is the number of neutrons scattered per unit volume 

4 ® = (9] 

and per unit time, but 





dE 2 E 



Thereforeo dE Z (E) . 

® dE s I 



» ; 

q dE 



E 



(15) 



E “ E ^ 
H 

E(l=a) 



is the probability a neutron will land in the X th groups 



JL'k 

y is the number of neutrons scattered in the k th group that land 

in the ^ th grot^Jo Since the absorption, sources and leakage in a 

large thermal reactor are small, the slowing down density q, is rather 

constant over an energy groupo This fact is used to compute the num= 

t “ [10] 

ber of neutrons scattered from the group k to group^ « 

The equations for each case are listed belowo 



Case 



lo E. > oE 
k 



, Equation 

s . _ w 1 1 

° (i"“> ' 



■) (16a) 



2e E - > o E 
-a k 



oE >E. >oE 
k JL+1 k+1 



jtk 



r 



(l<=>a) ln(E /E 

k k+1 



E.-E. 

/ ^+1 



E 



k|l 



jg 

4^ a « ““ 

E 






(16b) 
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(16c) 



3o 

E„ < a E 
J^+1 k+1 



^Lk 



r 






(1-a) 





1 


OC 


'E 




. ““ 


■ A 


E 


E 


L 


_ k+1 


k J 



E 



1 a In 



k+1 






E^ < a E 
Ji*l k+1 



Zk 



»y^ te ' ■"' ■ ' ■" 

^ (1-a) In 









(I6d) 



la 

“ — - ~ ” ® In — r“ 

E E. a E 

. k+1 J k+3j 



“ \*1 



JLk 



r * 0 



(I6e) 



For tho cases where the scattering occurs within the energy group 
we find 



E , > a E . 
/+! ^ 



E- < a E 
A*1 JL 



. A '2 r E - E 

JLl I 3L + ^ 



Ji-i 

y 



[- 



in 



(I6f) 



X 



^ ^ *1 



We may now proceed to evaluate equation i| by integrating over the 
r and z variables# However^ before doing so we must consider the direction 
in which integration should be carried out* It is always best to integrate 
in the direction of neutron flow since any errors produced in the process 
of integration will be diminished as we proceed further in this direction® 
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If integration is carried out in a direction opposite to that of neutron 

[7] 

flowp any errors produced will be magnified* 

We must now difference our equations with respect to r and z* 

We can arrange terms in equation (8) in order to facilitate the Integra^ 
tion* Thus equation (8) becomes 



[ 



\ ■^] ■ -[«f 





( 17 ) 



We must now further divide the results of the integration into two 
separate cases? 

Case A <=> The component of velocity in the r direction is less than 

the component of velocity in the z direction* 

Case B - The component of velocity in the r direction is greater 

than the component of velocity in the z direction* 

If we consider Case and difference r and z ais illustrated in the figure 

belowg we obtain four distinct subcases* These will be classified in 

accordance with the direction of 

travel being determined by the signs 

of b and h * A positive sign 
i f 

indicates a direction of travel away 
from the center of the core* The 
various cases are listed in Table I* 



A-9 




Tabl« I 



In all A cases |b^j < jh^ | 
In all B cases I bj >|h^| 



Subcase 


b 

i 


h 

f 


Al, B1 


- 


- 


A2, B2 


•f 


- 


A3s, B3 


+ 


♦ 


Alip Bii 


- 


♦ 



Considering subcase A-3 and integrating equation 5 and knowing 

that 



dr 

" 3-1 





Figure IV 

Position Mesh for Case A 




and 

®IB-1 

where A = r 
5 3 

we obtain the 



^ fj 

m 









- r 

j-l and 6 * z - z , 

111 in m-1 

following equation 



A»10 




Where 2?^ ^ is the total cross section of the material at the center of 
the grid in figure lij i.Co^ the material to the left and below the j«>m 



point and r 



d-1 



is j >1 



The radial direction of integration is determined by the sign of p 

ioSoj by the sign of the b^ or d^* Consequently, the signs of , b^, 

d,, and A all change together* Likewise, the direction of integration 
^ J 

axially is determined by the sign of p, i*e*, by the sign of h © 

I f 

Consequently the signs of p., h and A all change together* If we 

f m 

use absolute magnitude signs for these quantities all Case A equation# 
will be identical insofar as the signs of various terms are concerned© 
Similarly, all Case B equations will be identical* 
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f 




If we now add the fpUowii^ expfessipn to both sides of equation (18) 



•>1 



. V^'*' 



we obtain 

H i\l 
M ' 1^.11 



;i,in 



V*" - V'"’^ 

j,m-l 






f j“l r 

0-1-* 



,3,ri '^5-l,ra-l 



^f yi 



+ — 
r 






-_i£L 






0-1 



j-l,m j-l,m-l 



Wl 
.^IS> 



jjin 0>ni-l 



'i«f 



• s +5 

0*>m 0 



i-l,m-lj 









j-l,m-l j-l,m-l j,m 



0,m 



1 1 Ini ^'*’ ^ ^^1^ 

If this expression is multiplied by A cO where A 1 o3 ““i 

'0 O’* |hJ 



m 



we obtain 

I b J (1-40*^) 



j,m j,m 0>ni-l J,m-1 



A; 

CJO^ 



\^ \ / A C 



3} 2 j-1 r 



+ A A)'^' 



(18a) 



(18b) 
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Where 




And finally tranaposinj^ terns in equation (l8b) we obtain an expression 
for the vector flux at the point j,ni* 




By use of sinllAi* procedures for developing the equations for the 
three remaining A subcases we obtain the three equations on the 
following pages* 
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For* subcase A«l| 




(19b) 



For subcase A-1 





A similar procedure is used in deriving the equations for the 



Figure V 
Position Me^h 




B cases; except that as shorn 
in Figure $ the r and z com- 
ponents of the flux are different 

because b > h • The method 
i f 

of deriving the flux equation 
for case B-3 is outlined on 
the following page. The 
equations for the other B cases 
are listed also. 
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V'* - ^ 







1^1 
* 3-1 








’f d-1 




2 

N 




V'^’^ ♦ 

0,m 0,m-l 

S + S 
j,in j-l,m-l 




SL^-x,t 



(22a) 



To both sides of equation (22a) we add 



-i', V’" * V‘^’' - V’" 

1 i . V 3-l.»>-l 3-l,m-l 3," 

0 » 






0>m 



and define 



«0 



•if 



tl 'll 

A 1 b 



ml X 
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V obtain 
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- V 



'■t. 






j-l,m 



J-l,ia 







♦ ^ ^ 



f J-1 



d-1 



V'" *A 

3 



(22b) 



Prom this expression 

9- 



lb 1 (T- J"") fV'" - - V'" - ^ 

_l i> 60 > V "i-l,ro \i-l,m >l«m-l .i-l,m-l, 

j,m \^ \ / /I c \ 

* .i / J?m-1 i A 



1+A. 



^ 2*^ I V 




,m-l 

yi 



+ — 
r 



d-1 



(22c) 



where k\±s defined by equation (l8c) 



For case B-U 



, i» 1 (1- - V'" V'"*^ 



+ A, 



;3,m 






(23) 



where A is defined by equation (19b) 

a 

For case B-1 



r, |bl(1- '"'")rV’^ « 

y ipf J il L j-*-l,ra .i-Htm'fl 1 



2 pf 3-1 



(2U) 



where A^ is defined by equation (20b) 
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For case B-2 



j,D 



1 - V-^’" 1* 



I il a 3-1 r ^ 



( 25 ) 



where A^ is defined by equation (21b)« 

The vector fluxes defined by equations (l8d), (l9a), (2ba)^ (21a)^ 
(22c), (23), (2k) t (25), are valid for all values of (x with the excep- 
tions of the end points of our ^l space where p ■ +1* At these points 
trouble is e:q)erienced with singulariJtio^^^However, a modified equation 
may be derived by letting the terms containing T in all equations 

following equation (l) be replaced by zero, since for p ■ + 1 the; 

' 1 

term iiir ecjuation (l) vanishes. This special equation was not 
needed since values of p were chosen which avoided these singular 
points. 

The procedures described above result in n equations. One more 
equation is needed and is provided by setting n ■ -1 in equation (3)» 

The developAient, above may be used when the values of b, c, 0 , 0 

® - 1 , 

and Q used above are replaced by -1, 0, J^(f^*-l), 0, and 1, respectively. 
These values reduce equation (6) to that found for the special case. 

Thus equation <(6c) becomes with i ■ 0 



j,m 






1 ♦ 






2 
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where 
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COB 






Nomenclature 



di 



Values obtained by differencing and integration of the transport 
equation over ^ , and defined by Equations 7a to ?d of Appendix 
A. ^ 

Values depending upon the quantized values of , and defined 
by Equations 7e and 7f of Appendix A* 



a \init vector in the direction indicated by a subscript* 
i 

Fraction of fission neutrons which are produced in a particular 
energy group Ji • 

r The radial distance to a point . 

r The mean radius between two adjacent grid points* 

h 



z The axial distance. 

E Neutron energy from a point* See Figure 1, Appendix A* 

Values obtained by the use of -ttie Gauss quadrature approximation 
to sum the flux over p, space, defined by Equation 13, Appendix A* 

Values obtained by the use of the approximation to sum the 
flux over p, space, defined by Equation 12, Appendix A, 

oc Ratio of the minimum neutron energy after scattering to its 
energy before scattering* 
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a 

Q 






■» 

2 

2 

f 

Q. k 

y 



0 (r,z 









'V 



o> 

/s 

-TL 



Ch*id spacing of a mesh. 

Cosine of the angle 

Angle between the axis and the direction in which neutrons 
travel* See Figure I, Appendix A, 

Cosine of the angle 0* 

I 

The probability that a neutron which is scattered in the k 
group will land in the group. 

Number of neutrons produced per fission* 

Total macroscopic cross section* 

Fission macroscopic cross section* 



The probability per unit length of neul^n Jravel that a 
neutron is scattered in group k* 

, Neutron flux -vdiich is classified according to speed, v, 
v) 

colatitude and azimuthal direction cosines p and||^, and 
position r and z* 

Colatitude angle denoting neutron travel with respect to the 
radius* See Figure I, Appendix A* 

See defining equation Page A-12, Appendix A* 

See defining equation Page A-17, Appendix A* 

Direction of unit vector in which neutrons travel* 
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Superscripts and Subscripts 



i Index labeling quantized values * 

f Index labeling quantized values of 

j grid position index in the radial direction* 

m grid position index in the axial direction* 

neutron speed index* 

k neutron speed index of group in which neutron is scattered* 
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I 




Prior to the iise of this program the following preset parameters 
must be defined and/or calculated* 

pzl « number of grid points in a radial direction less one 

pz2 * number of grid points in an axial direction less one 

pz3 * number of iterations, N, allowed at a set radius, in 

attempting to achieve a convergence of source’ 'pat tenu--. 
pzli ■ number of adjustments, t.N'^ of radiua_Allowed in 
attempting to achieve criticality* 
pf2 ^ drum address of cross section 
pfl + 92(pzl ♦ 2) (pzl + 2) 
pf3 » 92(pz2 +2) 
pflt * 6(pz2 + l) 

The preset parameters listed below were used in the developanent of 
the program and have the values as indicated 

pfl ■ 5000; drutn( address of first flux at origin point* 
pf5 “ 70lij address in fast memory at which the print-out 
routine will be placed when needed 

V 

pflO “ 2650} address of the print-out routine in the drum 

pfl5 = 300j length of the print-out routine 

pfl3 “ lli70j address in the fast memory of the vector flux 

subroutine of the quadrant being calculated* 

It is necessary also to furnish the program at 

l) a9, the value of +(2pzl)*0* This value is used in the 

A ^ 

calculation of ^ 



2) rli, the value of the first trial radius of core 

3) r3> the value of the second trial radius of core 
For a three group solution, values of each of the thirty-six 

fluxes and the three sources at each point are assumed and placed on 
tape in the order indicated in Table I. 

Table I 

Sequence of fluxes at a point 



f i / 

0 0 0 

1 

2 

1 0 

1 

2 

2 0 

1 

2 

10 0 



0 .< f < 3 a <i < 2 0‘ kj( < 2 
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The points are ordered as indicated in Table H* 



Table II 

Sequence of points 



d 



0 



m 

0 

1 

2 



1 



pz2 + 1 
0 



pa2 * 1 



pzl+1 



pz2’fl 



As explained in the Procedure the probability of scattering 
from group to group and fission cross section are associated with the 
grid points. These are also placed on the tape of data immediately 
following the flux for f * 3# i * 2, and * 2 in the order 
indicated on the following page in Table HI* 



Table III 



Sequence of additional data for a point 



T 



11 



r 

r 

Y 

7 



21 



31 



22 



r 



32 

33 



Sf 

The total scattering cross section is associated with the region 
between points and are ordered as in Table ]V* 



Table 17 

Sequence of Total cross sections 

ji S 

0 0 

1 






pz2 



2 



pzl 



0 

pz2 



fc TAPE Iil0-316-1032 S AND GAUSS CAMPBELL AND PANCIERA 

(2U,6) 

pz2«5 pz3=20 pzl;*20 pf3,*5000 pf2«lOOOO pf3*6I|li 

pf5-70li pfl0«2650 pfl3-lit70 pfl?-300 



pzl-5 

pfli'»36 

/Constants 

al,+l*0 

a3,+2,li6 



-0*333333 +0*6666667 a2,+0*0 +1*333333 +2*0 

aJt,-l*0 a^,+l*9677051i +1.1306938 

+1.306938 +l*9677051t a6,+0*03 +O.03 +10.0 a7,+0^3l83 

+o;363it +0*3183 a8,+0*652lU^2 +0*3U7851|8 +0*3l»78^1i8 

+0*6^21ii52 a9,+22*0 al0,+U*0 aH,*l,69hh62$ +0*381*33la 

/Cross Sactions 



bl<»+l*0 DITTO 12r/\)2,+i;o 

m85pf3 m9,pfl* mlO,+92 mil, +6 

ml5,+0 
/RADII 

r3j,+0*0 rlt,+0i0 r$,*0*0 

/NORMALIZING CONSTANTS 



DITTO 12r/m6,+0 m7,+0 

ml2,p£l ml3,pf2 mllt,+0 



dl^+0.0 

/VARIABLES 


d2,+0*0 


d3,+0.0 rl,+0.0 


r2,+0*0 


sl0,+0*0 




Cl3+li*8 


c2,«5*6 


c3>+0.688^ 


+0.229^ 


+0*1;590 


+3.0355 


+1*0118 


+2.0236 










/j~l,m-l numbers 










nlj+l.O 


DITTO 


72r/sl,+l*0 


+1*0 


+1.0 


1PL,+1*0 


+lio 


+1.0 


+1.0 +1^0 


+1*0 


+1^0 




/3p*lj)m numbers 










n2,+l«0 


DITTO 


72r/s2,+l*0 . 


+1.0 


+1.0 


t2,+1.0 


+1*0 


+1.0 


+1*0 +1*0 


+1.0 


+1*0 




/j-=l,m+l, numbers 










n3,+1.0 


DITTO 


72r/s3,+l*0 


+1*0 


+1.0 


t3,+l*0 


+1.0 


+1.0 


+1*0 +1*0 


+1*0 


+1*0 




/o, m“l, numbers 










nit, +1.0 


DITTO 


72r/sli,+l*0 


+1*0 


+1*0 


tit, +1.0 


+1.0 


+1.0 


+1.0 +1*0 


+1*0 


+1.0 




/j,m numbers 










n^,+l*0 


DITTO 


72r/s$,+1.0 


+1*0 


+1*0 


t5,+i*0 


+1*0 


+1.0 


+1*0 +1*0 


+1*0 


+1*0 




/j<>m+l, numbers 










n6,+l*0 


DITTO 


72r/s6,+1.0 


+1*0 


+1.0 


t6,+l*0 


+1^0 


+1*0 


+1*0 +1*0 


+1*0 


+1.0 





/j+lpiS“l 5 numbers 



n7j.+loO 


DITTO 


72r/s7,+1.0 


♦1.0 


♦1.0 


t7,+1.0 


+loO 


+leO 


♦1.0 


♦;.o 


♦1.0 


♦1.0 




/j+lam^ numbers 












nSj+loO 


DITTO 


72r/s8,^1.0 


♦1.0 


♦1.0 


t8,+1.0 


♦1.0 


+1«0 


♦1.0 


♦1.0 


♦1.0 


♦1.0 




/j+lj,m+l, numbers 












n^j+l.O 
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